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Abstract 

A new version of the two-step multi-boson algorithm is developed with 



^ ■ different fermion actions in the multi-boson and noisy Metropolis steps. 



1 Introduction 



O ' Multi-boson algorithms for dynamical fermions |1|, 0, 0, ^ offer promising possibilities for 
I numerical simulations of QCD and other similar theories The two-step multi-boson 
^ [ (TSMB) algorithm has been succesfully applied, for instance, in super symmetric Yang- 
D '. Mills theories B and for simulations of QCD with SU(2) colour at high densities R]. In 
the present paper I shall consider a TSMB algorithm allowing for different actions in the 
' steps. 

^ ' For convenience of the reader let us first shortly summarize the main features of 

■ ■ ■ ' TSMB. Let us consider the case of an arbitrary number of identical fermion flavours Nf 
and assume the existence of a hermitean fermion matrix Q, which has the determinant 
det(Q) appearing in the effective action for the gauge field after the integration over 
the fermionic variables in the path integral. Multi-boson algorithms are based on the 
representation |l|] 

|de.«)r'^{det,g'^)}"-.-^. ,1) 
Here the polynomial P„ satisfies 

lim P„(x) = x~^f/^ (2) 
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in an interval x G [e, A] covering the spectrum of Q^. For the multi-boson representation 
of the determinant one uses the roots of the polynomial rj, {j = 1, . . . ,n) 



(3) 



Assuming that the roots occur in complex conjugate pairs, one can introduce the equiv- 
alent forms 

n n 

Pn{Q') = roUKQ ± + = ^0 n (0 - P*)iQ - P,) (4) 

i=i i=i 



where Tj = (fij+iuj)"^ and pj = fij+iuj. With the help of complex boson (pseudofermion) 



fields one can write 



N 



det(g) oc Yldet[{Q-p;){Q-p,)r' 



oc / m exp ^ - E E K HQ - P*,)iQ - P,)]y^ 

j=l xy 



(5) 



Since for a finite polynomial of order n the approximation in (|^) is not exact, one has to 
extrapolate the results to n ^ cxd. 

The difficulty for small fermion masses in large physical volumes is that the condition 
number A/e becomes very large (10^ — 10^) and very high orders n = O{10^) are needed 
for a good approximation. This requires large storage and the autocorrelation of the 
gauge configurations becomes very bad since the autocorrelation length is a fast increasing 
function of n. One can achieve substantial improvements on both these problems by 
introducing a two-step polynomial approximation: 



X 



?i2— >oo 



X G [e. A] . 



(6) 



The multi-boson representation is only used here for the first polynomial P^J-* which 
provides a first crude approximation and hence the order ui can remain relatively low. 
The correction factor P!^^ is realized in a stochastic noisy Metropolis correction step with 
a global accept-reject condition during the updating process. In order to obtain an exact 
algorithm one has to consider in this case the limit ^2 — oo. 

In the two-step approximation scheme for Nf flavours of fermions the absolute value 
of the determinant is represented as 



det(g) 



Nf 



1 



detPiP(g2)detPA2^(Q2) 
2 



>(2), 



(7) 



After an update sweep over the gauge field a global accept-reject step is introduced in 
order to reach the distribution of gauge field variables \U] corresponding to the right hand 
side of (0). This can be done stochastically by generating a random vector r] according 
to the normalized Gaussian distribution 

and accepting the change \U] [U'] with probability 

PaHU'] - [f/]) =min{l,A(r/;[f/'] ^ [U])} , (9) 

where 

MV, [U'] - [U]) = exp {-v^Pi'J{Q[Uf)r^ + Pi'J {Q[Ur)v} ■ (10) 

The Gaussian noise vector r] can be obtained from rj' distributed according to the 
simple Gaussian distribution 

J[drj']e-^"^' ^^^^ 

by setting it equal to 

V = PH\Q[Ur)-'^v' ■ (12) 

In order to obtain the inverse square root on the right hand side of (p!2D, one can 
proceed with a polynomial approximation 

PiJix)c:^P,^Jix)-K XE[0,\]. (13) 

This is a relatively easy approximation because Pj^\x)^^ is not singular at x = 0, in con- 
trast to the function x"^-^/^. A practical way to obtain P^^^ is to use some approximation 
scheme for the inverse square root. A simple possibility is to use a Newton iteration 

The second term on the right hand side can be evaluated by a polynomial approximation 
as for P(2) in (I) with Nf = and P(^) pf ^p(2)_ ^jj^g iteration in (|T|) is fast converg- 
ing and allows for an iterative procedure stopped by a prescribed precision. A starting 
polynomial Pq^^ can be obtained, for instance, from P^^^^ in (|^) with Nf —>■ —^Nf. 

The TSMB algorithm becomes exact only in the limit of infinitely high polynomial 
order: ^2 ^ oo in (|^). Instead of investigating the dependence of expectation values 
on n2 by performing several simulations, it is better to fix some relatively high order n2 
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for the simulation and perform another correction in the "measurement" of expectation 
values by still finer polynomials. This is done by reweighting the configurations. 

The reweighting for general Nf is based on a polynomial approximation P^^-* which 
satisfies 

Ji^^P!i^{x)Pil\x)Pif{x) = , X G [e', X] . (15) 

The interval [e', A] can be chosen by convenience, for instance, such that e' = 0, A = Xmax, 
where Xmax is an absolute upper bound of the eigenvalues of Q^. In practice, instead of 
e' = 0, it is more effective to take e' > and determine the eigenvalues below e' and 
the corresponding correction factors explicitly [^j. With properly chosen [e'. A] the limit 
71-4 — * oo is exact on an arbitrary set of gauge configurations. For the evaluation of P^^-* 
one can use ?T,4-independent recursive relations P], which can be stopped by observing the 
required precision of the result. After reweighting, the expectation value of a quantity O 
is given by 



{OexpW{l~p(f{Q^M)u,, 
{e^p{vni-Pifm]v})u,, 



(^) = —, r-iz — ^(4),^... ' (16) 



where is a simple Gaussian noise like rj' in (|TT|). Here (. . .)[7_^ denotes an expectation 
value on the gauge field sequence, which is obtained in the two-step process described 
above, and on a sequence of independent ?7's. 

The necessity and importance of reweighting the gauge field configurations depends 
on the condition number A/e, which grows roughly as the squared inverse of the fermion 
mass in lattice units. For relatively heavy fermions one can choose n2 high enough, such 
that the effect of reweighting is negligible compared to the statistical errors. This can be 
checked on a subsample of statistically independent configurations and then the systematic 
errors due to the finiteness of n2 are under control. The reweighting becomes necessary 
only for very light fermions as, for example, in [0]. In cases if reweighting becomes 
important the computational work for obtaining the reweighting factors is comparable for 
the calculation of the inverse of on the vectors i] in (|16D . This is typically less than the 
off-line calculations performed, for instance, for determining some correlators and their 
matrix elements. Of course, if the reweighting has some effect, then it has to be taken 
into account in the process of estimating statistical errors (see P|). 



2 Variable multi-boson update step 

The TSMB algorithm is based on the fact that the change of the fermion determinant 
in the update sequence can be approximated by a substantially lower order polynomial 
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than would be required for the final precision of the simulation. More generally, one can 
also allow for the fermion action in the multi-boson step to differ from the true action 
describing the model to be simulated. (This has been done in a special case in ref. ||T0| 



where the update step is performed with the pure gauge action.) Of course, the correction 
steps become in general more involved because one has to correct for the difference of the 
two actions, too. 

Let us denote the auxiliary hermitean fermion action in the multi-boson step in general 
by Qo. Instead of (^ the fermion determinant is now represented by 

~ 1 

det(Q) ~ wrp-^ -...x TTy--^ — . (17) 

detPffiQl) det P£)(Q§) det pZ\Q') 

The main difference compared to (|^) is that here an additional polynomial (P*^^^) appears 
which is needed in order to compensate for the use of a different action Qq in the first 
polynomial P^^\ The polynomials appearing in ([T7|) have to satisfy, with R^"^ ~ 1 and 
R^"^ ~ 1, the relations 

^-""^ - Pi^.'m.) - PV'H.) ■ ^'"wp<'^'w^fi-'w. (18) 

Here i?*^^-*(x) ^ 1 and R^^'^\x) ~ 1 have to be much better approximations than R^-^\x) ~ 
1. 

The factor (det P*^^^)^^ in ([TtD is generated by the multi-boson update step using the 
action Qq. The remaining factor, which has to be reproduced by the noisy Metropolis 
step, can be written, for instance, as 

1 1 

^ = = — = fl9) 

det P(2)(g2) detP(2)(g2) det(p(2)(Q2)i detP(2)(g2) p(2)(g2)i^ • 

Note that here and in what follows the roles of P*^^^ and P*^^-* can also be exchanged. 

Using the form in ([19D one can write the detailed balance condition for the acceptance 
probability Pa in (H) as 



^ lu]) _ det (p(2)[[/]lp(2)[t/]p(2)[[/]| 



Pa {[U] ^ [U']) det (p(2)[f/']^p(2)[t/']p{2)[[//]f 

/ [d?]] exp (-r/tp(2) [t/'] 5P(2) [[/']p{2) [Uf2n^ 
/[dr/]exp(-r7tp(2)[f/]|p(2)[f/]p{2)[f/]i^) ' 

Here the gauge field dependence of the polynomials is displayed, which refiects the de- 
pendence of the fermion actions on the gauge field. The detailed balance condition can 



be satisfied similarly to (P) with 

A{r]; [U'] ^ [U]) = exp [-7]^ P^^^[U']^ P^^^[U']P^^^[U']^7] + 7]^ P^^^U]'^ P^^\U]P^^^[U]'^r]} . 

(21) 

The required Gaussian can now be obtained from a simple one in ([TI|) by 

r^ = p(2)[f/]-|p(2)[[/]-^' . (22) 

For the evaluation of ( pT]) we need the polynomial approximations 

Pif(x) ^ Pifix)-'^ , Pifix) ^ Pi?(x)-^ , X G [0, A] , (23) 

and 

P(-33)(X) ^ Pif (x)-^ ^ Pi^)(x)^ , Pt_!\x) ^ Pff (x)-i Pi?(x)^ , X G [0, A] . 

(24) 

The measurement correction can be performed in an analogous manner as before. For 
this we need the polynomial approximations P^^\x) and P^^\x) corresponding to ([T8|) : 

- Pi^>mm) • - . (25) 

The final precision is given now by the quality of the approximations R^'^^\x) ~ 1 and 

^(124) ~ 1. 

3 An application 

Apart from the doubling of the number of polynomials in the noisy Metropolis correction 
step the present algorithm is entirely analogous to normal TSMB. As a first application 
let us here consider the case where the two hermitean actions Qo and Q are the same 
Wilson fermion action, except for the values of the parameters /? (the gauge coupling) 
and n (the hopping parameter). 

Let us denote the parameters in the multi-boson step by (/3o, /^o) and in the noisy 
Metropolis step by (/3, k). An interesting question is how for fixed parameters (/5, k) the 
change of (/3o, ^o) does influence the behaviour of the algorithm. To see this I performed 
a test run with two flavours of Wilson fermions on an 8^ ■ 16 lattice at parameters [jS = 



5.28, K = 0.160). (This corresponds to a point on the Nf = A cross-over line of fill.) The 
main polynomial parameters were as follows: [e. A] = [0.00875, 2.8], rii = 24, n2 = 70. 

The change of the acceptance rate in the noisy Metropolis step for changing /3o resp. 
kq is shown by figure |l]. As one can see, the acceptance remains quite good in a relatively 



broad range of parameters. This is partly due to the fact that the distribution of the 
exponent in the noisy Metropohs step becomes substantially broader when the multi- 
boson update parameters {Pq, kq) differ from k) (see figure 0). 

Changing (/3o, kq) can be used to introduce more randomness in the update process. In 
fact, since the final distribution of gauge configurations does not depend on the parameters 
of the multi-boson update, during a run one can randomly change (/3o, k,q)- This improves 
the autocorrelations. As an example, the measured integrated autocorrelation of the 
plaquette in the original TSMB run with the above polynomial parameters is rfj^^'^ = 
2.9(3) ■ lO^MVM (Matrix- Vector Multiplications by the even-odd preconditioned Wilson- 
Dirac matrix). This is decreased to rfj^'^'^ = 1.2(3) ■ lO^MVM if (for fixed kq = k) (3q is 
changed randomly during the gauge update in the interval (3q G [5.16, 5.40]. For counting 
the number of MVM's in an update cycle the following approximate formula is used: 

Nmvm ^ Gin^Ns + Ng) + 2Nc{n2 + ^3 + 713 + n_^) . (26) 

Here Nb is the number of boson field updates, Nq the gauge field updates and Nc the 
number of noisy Metropolis accept-reject steps. The autocorrelations were determined 
in independent runs which were at least as long as hundred times the measured auto- 
correlation. In these test runs with relatively heavy quarks the polynomial order n2 is 
high enough and the effect of reweighting in (0) is much smaller than the statistical 
errors. Therefore, the autocorrelations were determined without taking into account the 
reweighting factors. 

The decrease of the plaquette autocorrelation as an effect of updating with a "wrong" 
action is at the first sight against intuition. This unexpected effect can be understood as 
due to partly compensating the dominance of the bosonic contributions against the pure 
gauge contributions in the effective gauge action. As a consequence of the gauge field 
part in (^Tj), the fluctuation of the gauge fields within the natural band of fluctuations is 
intensified. At the same time there is some decrease of the acceptance, too, but in total 
the effect of the amplified fluctuations is more important. 

An interesting question is the volume dependence of the effect of action variations 
in the update. In general, the allowed action changes may be expected to decrease in 
larger lattice volumes. However, one has to observe that the range with good acceptance 
in figure |I| is much larger than the range which would be allowed for entirely updating 
with different parameters and performing the necessary reweighting afterwards on the 
equilibrium configurations. The main reason is that the distance between two noisy 
Metropolis steps is typically less than 1% of the plaquette autocorrelation distance. In 
addition, on larger volumes the number of boson fields rii is also larger and the dominance 
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of the bosonic part in the effective gauge action is stronger. Therefore the question of the 
volume dependence is not easy to answer and can be best done by performing actual test 
runs on larger volumes [1^ . 



Another interesting possibility is to apply TSMB for variable actions in numerical 
simulations with improved fermion actions. Using a simplified version of the particular 
improved action in the multi-boson update step facilitates the implementation of TSMB 
algorithms and may also improve the autocorrelations. 
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Figure 1: The average acceptance of the noisy Metropohs step at (/3 = 5.28, k — 
0.16). On left (right): changing Pq (kq) in the multi-boson update step. 
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Figure 2: The distribution of the exponent in the noisy Metropohs step at {(3 = 
5.28, K, = 0.16). The parameters in the multi-boson update step: on the left (right) 
Po = 5.28, Ko = 0.160 {Po = 5.28, kq = 0.157). The vertical lines show the mean 
value of the distributions (and zero). 



10 



